doseresNMA antidep/drnma.model.metareg.R

In htx-r/doseresNMA:

drnma.metareg <- function(){
  for(i in 1:ns){ # loop through ns trials
    w[i,1] <-0
    delta[i,1] <- 0
    u[i] ~ dnorm(0,0.001) # baseline effect

    for (k in 1:na[i]){ # loop through na arms within a trial
      r[i,k] ~ dbin(p[i,k], n[i,k])
      logit(p[i,k]) <- theta[i,k]
      theta[i,k] <- u[i] + delta[i,k]+(gamma[t[i,k]] - gamma[t[i,1]])*cov[i]

      # arm deviance
      rhat[i,k] <- p[i,k] * n[i,k] # expected number of events
      resdev[i,k] <- 2 * (r[i,k] * (log(r[i,k]) - log(rhat[i,k])) + (n[i,k] - r[i,k]) * (log(n[i,k] - r[i,k]) - log(n[i,k] - rhat[i,k])))
    }

    # study deviance
    resstudydev[i] <- sum(resdev[i, 1:na[i]])

    for(k in 2:na[i]){ # treatment effects
      # RE delta's
      delta[i,k] ~dnorm(md[i,k], prec[i,k])
      DE[i,k] <- ((beta1[t[i,k]] * dose1[i,k]) + (beta2[t[i,k]] * dose2[i,k])) - ((beta1[t[i,1]] * dose1[i,1]) + (beta2[t[i,1]] * dose2[i,1]))

      # multi-arm corrections
      md[i,k] <- DE[i,k] + sw[i,k]
      prec[i,k] <- prec.d *2*(k-1)/k
      w[i,k] <- delta[i,k] - DE[i,k]
      sw[i,k] <-sum(w[i,1:(k-1)])/(k-1)
          }
  }

  #** PRIORS

  # to delta's RE model
  tau.d~dunif(0,5)      # common standard deviation is given a vague prio
  prec.d<-pow(tau.d,-2) # common precision in RE-NMA model

  # to beta's
  beta2[ref.index] <- 0
  beta1[ref.index] <- 0
  for (k in c(1:(ref.index - 1), (ref.index + 1):ndrugs)){
    beta2[k]~dnorm(0, 0.001)
    beta1[k]~dnorm(0, 0.001)
  }

  # to gamma RE model
  gamma[ref.index] <- 0
  for (k in c(1:(ref.index-1),(ref.index+1):ndrugs)){
    gamma[k] ~ dnorm(g,prec.g)
  }
  g~dnorm(0,0.001)
  tau.g~dunif(0,5)      # vague prior to gamma's SE
  prec.g<-pow(tau.g,-2)

  # overall deviance
  totresdev <- sum(resstudydev[])


  #** drug absoulte effect
  # baseline effect
  for (m in 1:ns){
    rr[m,1] ~ dbinom(p0[m],nn[m,1])
    logit(p0[m]) <- OO[m]
    OO[m] ~ dnorm(Z, prec.z)
  }
  # baseline priors
  Z ~ dnorm(0, 0.001)
  prec.z <- pow(sigma.z,-2)
  sigma.z ~ dunif(0,5)
  p.placebo <- exp(Z)/(1+exp(Z))
  # combine baseline with OR -> drug absolute effect
  for (k in c(1:(ref.index-1),(ref.index+1):ndrugs)){
    for( j in 1:nd.new){
      OR[j,k] <- exp(beta1[k]*new.dose[k,j]+beta2[k]*f.new.dose[k,j]+g*cov.pred)
      odds.drug[j,k] <- OR[j,k]*exp(Z)
      p.drug[j,k] <- odds.drug[j,k]/(1+odds.drug[j,k])
    }
  }


}